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Large eddy approximation of turbulent flow in DSC schemes 


Abstract. Large eddy simulation of turbulent flow is given a natural setting within the DSC framework 
of computational fluid dynamics. Periodic cellular coarse-graining prevents the nodal flow from piling up and 
preserves its large patterns. The coarsening operations are consistent with the near-field interaction principle of 
the DSC scheme and are just as simple in the treatment of boundary conditions. Numerical examples validate 
the approach. 
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L’acte de la nature est simple, et la complexité qu’elle parait avoir construite 
n’est que l’entrecroisement sans fin des antagonismes qui se sont neutralisés 


Henri Bergson 


1. Introduction 


Discrete schemes have to cope with the cursed situation that they first create by breaking space into 
pieces. Mesh cell systems, however fine, can never be perfect substitute for continuous space. Already 
in linear algorithms they produce artefacts, such as unphysical spurious solutions, and more harmful 
things happen in the non-linear case. 

In fluid dynamics, when the non-linear spectral transfer properties of the Navier-Stokes equations 
(the energy cascade; cf. PopE [Po]) come into play beyond the transition to turbulence, eddies are 
locally excited, down to very short scale; cf. KoLMOGOROV [Ko]. Eddies smaller than the cell size 
cannot be properly resolved by a mesh of realistic coarseness and thus tend to induce local fluctua- 
tions that artificially pile up. Large eddy simulation (LES) aims to discard such divergences through 
regularizing the flow by means of a suitable averaging filter, and to screen its essential large patterns 
in this way. 

LES stands for the most promising line in modern turbulence modelling. The name has been coined 
by DEARDORFF in 1970, who first applied LES methods to turbulent channel flow [De]. Important 
elements had yet been kept ready, formerly. Already REYNOLDS used temporal flow averages that 
as Reynolds averages play still a rdle in conventional turbulence models. In 1922, RICHARDSON 
proposed mesh cell averages for smoothing down local fluctuations. Such so called box filters stand 
today for one line of LES. Spatial filters of varied type indeed classify modern LES methods. Excellent 
insight into the state of the art is gained from the recent book of BERSELLI et al. [BIL], which not 
only exposes the elaborate mathematical framework but also addresses many items of conventional 
turbulence modelling (the prominent /— € model of LAUNDER and SPALDING [LS], e.g.). 

The present paper differs from most related work in the field in that it sharply separates the natural 
and technical aspects of turbulent instability and strictly focusses on the latter within the computa- 
tional context in hand. The natural part, viz. essentially (modulo discretization) the Navier-Stokes 
equations for viscous incompressible flow, are thereby taken as a physical given. This is in obvious con- 
trast to widespread use: In many conventional turbulence models, as in some lines of LES, turbulent 
pile-up is controlled by modifying the underlying dynamic equations or even the physical constants. 
Such models use, for instance, Reynolds number dependent global parameters - usually an eddy viscos- 
ity, e.g. - or time and space filtered Navier-Stokes equations or any (e.g. stochastic) variant. All this is 
allowed, of course, as long as it provides numerical data in harmony with observation. Also, a profound 
change of the inner dynamics can very efficiently delimit turbulence and thus lead to a convergent 
scheme. It is yet, obviously, tuning the underlying physical phenomenon rather than eliminating the 
technical cause of divergences that are manifestly an artefact of discretization. (Turbulence is still 
enigmatic in many respects, but divergence is certainly not a feature of turbulent flow - observed in 
the continuum of nature, e.g.) Therefore, the central point of attack in discarding turbulent pile-up 
will be in this paper the weaknesses of discretization, while the Navier-Stokes equations are directly 
discretized as taken from physics - which not only combats the desaster at the origin, but also yields 
validated results. 
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2. DSC schemes 


A mesh cell system is not just poor surrogate continuous space, but it always introduces its proper 
structure. Cellular meshes, with that we are dealing here, artificially impose local cell-boundary duality 
upon space - and it is this simple fact that dual scattering channel (DSC) schemes match in a 
canonical fashion. 

DSC schemes are finite volume methods of a very specific type. In technical detail, they are exten- 
sively dealt with in [Hel], so it is sufficient here to recapitulate their essential features in a crash course 
like manner. The algorithm is characterized by a two-step cycle of iteration which alternately updates 
the computed fields within cells and on their interfaces. If the updating instructions are explicit, then 
a near-field interaction principle gives rise to a scattering process interpretation known as Johns’ cy- 
cle [Hel]. A pair of vectors that represent the same field within a cell and on its surface essentially 
constitutes a scattering channel. Equivalently, scattering channels are sometimes defined as pairs of 
distributions that ’measure’ the field within a cell and on one of its faces. The well known primal 
DSC scheme is JOHNS’ TLM algorithm, wherein scattering channels are visualised as transmission 
lines [JoB]. In many applications, such as computational fluid dynamics, the TLM picture of wave 
propagation fails, but one can still go back to the far more general DSC setup. 

Within any physical interpretation, cf. , each scattering channel represents a pair of scalar 
or vector valued distributions z? = (p, Z) and z” = (p™~, Z) which evaluate a physical field Z on 
the surface and in the interior, respectively, of a mesh cell. A distribution p pertinent to a cell face is 
called a port and p™ its nodal image, and the two are related by pull-back, viz. 


(1) (p42) =] pee 7) = (p72 oe); 


for every Z (of class C™, e.g.) , where o denotes the spatial translation o : R? — R® that shifts 
the geometric node (i.e. the cell centre) in the (centre of) the respective face. 

It follows that there exists a scattering channel representation of DSC states: If M denotes the 
mesh cell system, O¢ the boundary of cell ¢ € M (0¢ is naturally identified with the set of ports 
with nodal image in ¢), and L? the span of the z@ (ie. the linear range space), then each state in 
the mesh permits a unique representation in the space 


se Pp p~ 
(2) P := Le LD eo LEDs 


with canonical projections 7?’" : P + P?" into the port and node components of cell ¢ . Further- 
more, there is a natural involutary isomorphism nb: P > P, 


(3) nb: (z?, 2? )1H (2? , 2”), 


called the node-boundary map, which hence maps P? onto P” and vice versa. (The cell index ¢ is 
omitted here and in the following without danger of confusion.) 

It stands to reason that the nodal images of different ports (e.g. on different faces of the same cell) 
can represent the same physical field in a node and may even coincide there as distributions - just as 
two ports pertinent to neighbouring cells obviously represent the same field on a common face, if they 
connect two channels on it. In fact, the scattering channel representation of DSC states in is in 
general highly redundant - which may be utilized for process parallelization [Hel]. 

Within the algorithm, the port and node components are updated at even and odd integer multiples, 
respectively, of half a timestep 7 and are usually constantly continued as step functions over the 
subsequent time intervals of length 7. 

Note that existence, not necessarily explicit construction or application, of a scattering channel 
representation characterizes DSC schemes, and that less redundant representations are actually used 
in many implementations. In fact, the scattering channel representation is basically a theoretical means 
for describing (and deriving) certain DSC properties - such as the following. 

A fundamental principle - closely related to the COURANT-LEVI stability criterion - is near-field 
interaction. It requires that every updated state of a node or face depends only on states (along with 
their history) in scattering channels connected to the respective node or face - the latter here being 
for once identified with its adjacent face, if such exists in any neighbouring cell. 

As a consequence of near-field interaction, every DSC process allows for an interpretation as a 
multiple scattering process in the following sense. 
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Let for any process z = (z?,2")(t) incident and outgoing fields z?, and z”,, be recursively 
defined as processes in P? and P”, respectively, by setting 2? (t) := 22y4(t-43) := 0, for 


t <0,and for 0 <t=mrT;meEN: 
ee 


5) 


zme(t+ 5) i= 2"(t+ 5) — nboa?, (t). 

Then, at every instant holds z?(t) = nboz?,,(t — $) + z},(t) and 

z"(t+ 4) = nboz? (t) + z%,(#+ $) . Also, near-field interaction implies that every state is 
only a function of states incident (up to present time t) on scattering channels connected to the 
respective node or face. 

More precisely, by induction holds: 


a? (t) := 2° (t) —ndboz?.. Ct - 


(4) 


Theorem. There exists a pair of functions R and C, defined on back in time running sequences of 
incident and outgoing fields, respectively, such that for every cell ¢ € M the process z? = mZ0z 
complies with 2 

(5) ze (t+ 5) = R( (Zin (t — HT) ) pear; wen) 


and the port process z? = 72 © z satisfies 


a Ze (E+ 7) = C( (ate (t+ 5 —HT))ntacs nen) - 


( ’ |’ short-hand for ’in any of the (7 or 2) cells adjacent to’ ) 
Remarks 


(i) The statements immediately imply that 20 out and ze in are themselves functions of states incident 


on connected scattering channels, since 


elou(t + 5) = RU(2% (t —HT))peacs nen) — nbo 22 ,(t) and 


n 
2b in (t) a a Erne: _ 3 7 HT) ) nip; nen) Se 
~ Tr 
— nbo 2 out (t = Cy 


(ii) R and © are called, respectively, the reflection and connection maps of the DSC algorithm, and 
a field excitation may be implicit in @, cf. [Hel]. 


(iii) Near field interaction implies computational stability, if the reflection and connection maps are a- 
passive , i.e. contractive in this sense [He2]. 


3. Coarse-graining 


Mesh cell systems import artificial defects into space - which yet sometimes carry their remedy in 
themselves, at least in part: For instance, the cellular mesh, by separating the virtually unresolved cell 
interior from the coarse-grained cell-boundary skeleton grid, necessarily fixes - and hence delimits - the 
scale of local resolution. A clever design always takes into account the coarsening effect of the cellular 
mesh - which, on the other hand, provides quasi free of charge a basis for large eddy approximation 
in the DSC setup of computational fluid dynamics. In fact, already intuitively DSC schemes go along 
with kind of natural large pattern approximation that - in a sense now made precise - associates cell 
averages to some nodal quantities. 

The following definition is independent of any particular application and therefore given without 
reference to fluid dynamics. 


Definition (coarse-graining) 
For any port p € 0¢ with nodal image p~ let B? denote the set of all ports on O¢, the nodal images 


of which coincide with p~ as distributions. Also, let w, € [0,1]; re B?, be a set of weights such 


that So epp Ur = 1. 
Then, every DSC field which in the scattering channel representation on each component q € BP C O¢ 
equals (z%, 2% ) with 


(7) 2 = ene Ure” 


is named an (in ¢ and with weights w, ) coarse-grained field (pertinent to p - or also to its nodal 
image p™ ). 
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The nodal state of every such field, coarse-grained in any cell, is hence a convex superposition - i.e. a 
weighted mean - of states that represent the field on the cell surface. Also, substituting a DSC field 
with any coarse-grained field that (obviously) leaves the cell boundary states unchanged is then called 
coarsening this field (in any component). 

Should such an innocent looking procedure untie the Gordian knot of turbulence ? Remember that 
our modest aim is repairing some defects of discretization - rather than resolving the fundamental 
questions of turbulence in the way mathematical physics should do that. What matters here is that 
in the DSC setup of CFD well-timed periodic coarse-graining of the flow efficiently prevents turbulent 
pile-up, all with retaining the large eddies. Note that the coarsening period should be taken as small 
as necessary to ensure stability and precision, but large enough against the time step (by one order of 
magnitude at least - as a rule of thumb): Clearly, coarsening inevitably also interferes with the flow 
dynamics. However - since the skeleton field on the cell boundaries is left unchanged - the perturbation 
is minimal, and in fact negligible, for suitable periods. Through removing only the turbulent fluctua- 
tions of very short range (below a scale given by cell size) which tend to artificially pile-up, coarsening 
regularises the flow to just such a mild degree that its essential large patterns are preserved. 


4. Viscous Boussinesq-incompressible flow 


The port and node distributions of a DSC algorithm can be finite integrals - as is the case for the 
TLM method, where finite path integrals over electric and magnetic fields are evaluated in a dis- 
crete approximation to Maxwell’s integral equations [He4-6]. In the case at hand, they are simply 
Dirac measures that pointwise evaluate the fields within the cells and on their surfaces, and compos- 
ites of Dirac measures which approximate the gradients of these fields; cf. section (5). 

In the Oberbeck-Boussinesq approximation all material fluid properties are assumed as 
constant - except fluid density which only in the gravitational term varies linearly with temperature. 
The energy equation of a Boussinesq-incompressible fluid of velocity w with thermal diffusivity a, 
heat source(s) q, and negligible viscous dissipation and heat capacity is then the convection-diffusion 
equation for the temperature T’, cf. [GDN], 

OT 
The Navier-Stokes momentum equations for (the k-th component of) a Newtonian fluid of dynamic 
viscosity 1, under pressure p, and in a gravitational field of acceleration g take the form [e.g. Ptk] 
Oy Oe ee ee ss 4 
(9) gz (oun) + div(oted) + gradpe = nAde + Oger. 
With 0. =const and 0(T)=0~6(T(t,#) — Tx), wherein 8B = 0 '0o/ OT, this becomes in 
the OB-approximation 


au d 
dig) 42 Ae 4 2) a Te ee. 


oo Ot Qo0 Qc0 


Integrating these equations over cell ¢ with boundary 0¢ and applying Gauss’ theorem to integrals 
over A = divgrad and div( ft) yields with time increment 7 the following updating instructions 
for nodal T and tv; , these quantities and q being averaged over the cell volume V¢ 


T(t+=) = 4 


(11) : ; 

So (agradT — Tt): dS 
Ve Jac 

and edn. 

tx (t+ 5) i= de + 7(B(T — To) Ge — ——— ) + 
Qc 

(12) se 

+o f (— gradti, — tpt): dS. 
Vo ac Qc 


In equations (Ii) and (12) the former updates (at time t — 7/2) of the nodal quantities enter the 
right-hand sides in the first line, and the last former updates (at time t) of the cell face quantities 
enter the second line. 
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The nodal values of T and tw, are updated at the reflection step - after coarsening the velocity 
field, e.g. with weights proportional to the pertinent cell face areas, whenever the coarsening period is 
attained - while the port quantities that enter the surface integrals at the right-hand sides are updated 
on the connection step of the iteration cycle. The next section outlines how we can proceed with that 
in a non-orthogonal hexahedral mesh. 


5. The non-orthogonal hexahedral cell 


The physical interpretation of a DSC algorithm associates a smoothly varying (e.g. in time and space 
C-) scalar or vector field Z to port and node states z? and z” of a mesh cell system. 

Let any hexahedral cell be given by its eight vertices. Define then edge vectors (,e)v=o,...,11, node 
vectors (,b).=0,1,2, and face vectors (,f),=0,...,5, using the labelling scheme of figure [Ia 


1 3 
po i= Z ed eee pw = 0,1,2 
13 (-1)* 
(13) and ,f := 7 (ae em re tr. 
A Laane Tesi20) e) t6= 0,...,5, 


with all indices understood cyclic modulo 12 and A denoting the cross product in R?. 


Fig. 1. Non-orthogonal hexahedral mesh cell. (a) Edge vectors. (b) Node vectors. (c) Face vectors. 


At every cell face v € {0,...,5} and for any given 7 € R, the following time shifted finite differences 
of Z in directions ,,b (js = 0,1,2) form a vector valued function 


2(-1)'( 2" 4-7/2 —- 14? y4) if = [e/2] 
( 2y41Z? — ap Z? ) tr if p A [o/2] 


([a] denotes the integer part of « € R). The time increments are chosen conform with the updating 
conventions of DSC schemes (as will be seen in a moment) and are consistent. In fact, in the first 
order of the time increment 7 and of the linear cell extension, the vector .V?Z in the centre point of 
face & approximates the scalar products of the node vectors with the gradient VZ. More precisely, 
let for a fixed centre point on face » and e € Ry; the e-scaled cell have edge vectors ,e~ := €,e. Let 
also ,V®" Z,, denote function for the e-scaled cell (with node vectors ,b~ = €,,b). Then at the 
fixed point holds 


(14) NW Et s= 


1 ~ 
15 b, grad(Z) > = ,b-VZ = lim lim —V? Z 
(15) < yb, grad(Z) > uo V lim lim =. as 
as immediately follows from the required C!-smoothness of the field Z. 
To recover the gradient VZ from in the same order of approximation, observe that for every 
orthonormal basis (,w))=0,....m-1 of R™or C™, and for any basis (,,b),=0,....m—1 With coordinate 
matrix BY = < ,u, ,b >, the scalar products of every vector a with ,,b equal 


m-1 _ 
(16) < ,b,a> =v <ypb,pu>< ,u,a> = Bi ow 
=:a8 (Gy=ery ey 


Note that we observe EINSTEIN’s summation convention at the right side of equation (16) and hence- 
forth - but without ever summing over indices appearing as left-hand subscripts! 
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From equation (16) follows 


(17) a, = Wop, with — (of) = ((en)*) 

This applied to the node vector basis ,,b and (15) yields the approximate gradient of Z at face v 
(18) Nag. = aN 2 

The scalar product of the gradient with face vector ,f” = < ,f, pu>, v € {0,1,2}, is thus 

(19) Sighs = ae 2 a ev os 


Continuity of the gradient at cell interfaces yields linear updating equations for Z? on the two adjacent 
faces. In fact, for any two neighbouring cells ¢, x with common face, labelled v in cell ¢ and & in x, 
continuity requires 


(20) eS 8, 


Substituting (19) for °S and *S and observing the time shifts in (4) provides the updating relations 
for Z? at the cell interfaces. To make these explicit, we first introduce the quantities ,z}" 

(c = 0,...,5; pw = 0,1,2) 

2(-1)' 2"), if pw = [e/2I 


(anti Z? — a2” )ie_r/2 else 


(21) y(t) r= 


which in virtue of (I) yields ,z? = (p, Z) = (p~, Zo ,07") = 


= 2 |Z .o', where ,o : n +4 p denotes the nodal shift pertinent to face 1. In particular 
(22) 2/9] (t) — 2 (—1)’ 12? |t ; 


which together with (27) is consistent for py £ [z/2] with 


1 
(23) ei (ttr/2) =~ 5 (ame + au2P)(t)- 
From ({14] 19} [21] 22]) follows that 
(24) wo | t+r = ist (2. |t+7/2 — 2 (—1)'olt/?! pa? (eee) 


n u/2 
= ish (2 |t+7/2 — 5lt/ 2h |t+r) : 


Continuity of Z, ie. (Z? = *Z?, with QO) PI) then implies 
: gh S2M(t47/2) + Xs¥ *2(t4+7/2) 
(25) Maer). = =——= a 
78 + (-1)ete Xs [6/2] 

For completeness we agree upon setting ee (t+T):= oe (t+7/2) for uw & [v/2] (although this 
contains a slight inconsistency in that continuity might be infringed; this can easily be remedied by 
taking the arithmetic mean of the two adjacent values). - In fact, our agreement doesn’t do harm, 
since any discontinuity disappears with mesh refinement. 

We dispose, hence, of a complete set of recurrence relations for z? (given z” by the former 
reflection step) which at the same time determine the field components on face v and their gradients 


(26) V4Zy = yD He 
Essentially this constitutes the connection step of the algorithm. 
Nodal gradients are similarly (and even more simply) derived, using 


. 
VIZ (tts) 2= (auiZ? — yZ?)(t); w= 0,1,2 
in the place of (14), and then again ([8). With the node and cell-boundary values and gradients of 
T and w the nodal updating relations for the latter are immediately extracted from equations (TJ 
(12]) in section [4]. For equation this is essentially (up to the convective term) carried out in [Hel], 
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section 5, and the procedure remains straightforward in the case at hand. Note that a well-timed 
LES coarsening routine, cf. section B] should be periodically carried out before the nodal step of the 
iteration cycle in order to obviate instablities called forth by the energy cascade [Po]. 

The updating relations thus obtained are explicit and consistent with near-field interaction (only 
adjacent quantities enter ). So, they can optionally be transformed into scattering relations for incident 
and reflected quantities (6) along the guidelines of sectionB]- with established advantages for stability 
estimates [He2]. 


6. Pressure 


Unlike compressible flow - wherein a thermodynamic state equation relates pressure to density and 
temperature - Boussinesq-incompressible flow is conserved by pressure acting like a potential against 
violations of the local flow balance. Hence, the pressure is directly coupled to the flow divergence by 
Poisson’s equation: 


(27) Ap = (@x/T) divi. 


The solution p for fluid velocities updated in the connection step provides a pressure gradient, which 
in the next reflection step - via equations - repairs sporadic violations of the flow balance. 
Systematic test computations have shown that additional divergence clearing, as proposed by the 
author in earlier papers, e.g. [He3], is not necessary (nor even profitable) - and may be rather on the 
debit of computational performance and stability. 

In integral form, after using Gauss’ Theorem, equation (27) becomes 


(28) rf gradp- dF = ox | ad-dF, 
aC a 


which can be solved by a Gauss-Seidel routine or by successive overrelaxation, carried out after up- 
dating the face velocities in the connection step. 

In each iteration, firstly the discrete cell boundary integral Ig¢ :=_f, ae u-dF is computed and 
explicitely - by solving a linear equation in Jg¢ and the cellular pressures - the unique p” that exactly 
compensates Ig¢ so that is satisfied in every cell. In a second run over the mesh, continuity of 
the pressure gradient is re-established at the cell interfaces by updating the face pressures along the 
lines of the last section (viz. taking Z as p there). The loop is reiterated until (usually after a few 
iterations) equations (28) hold with sufficient precision. 


7. Conclusion and Completions 


The DSC approach to computational fluid dynamics allows for simulation of turbulent flow in a plain 
and natural way. Cellular coarse-graining efficiently prevents turbulent pile up and preserves the large 
eddies. 

Since the energy cascade starts long before transition to turbulence (actually it is effective over 
the entire energy spectrum [Pol, pp. 72 ff.] ), a well-timed coarsening routine should escort every 
DSC fluid flow algorithm. Periodic coarse-graining significantly improves algorithm stability, even in 
the laminar regime - this without overly altering the flow dynamics, if the coarsening period is chosen 
as outlined in sectionB] Handy criteria for a good choice - in relation to mesh refinement, e.g. - should 
be subject to further study. 

Details of implementation and special applications are clearly beyond the scope of the present 
study. However, we leave off with some examples for the purpose of illustration. 

The graphics displayed have been computed with author’s test program DANSE, Hern {He7], 
which combines a Maxwell field TLM Solver with a DSC heat transfer and flow algorithm based on 
this paper. 

The computations for coaxial RF power transmission line RL100-230 (with data realistically taken 
from an ion cyclotron resonance heating experiment in plasma physics) are in excellent keeping with 
reference computations (FLUENT) and empirical measurements. For a line with inner conductor made 
of copper and air dielectric at atmospheric pressure the heating process has been simulated in horizontal 
position from standby to steady state operation with 160 KW transmitted CW power CW power at 
100 MHz frequency, with outer conductor cooled at 40 degrees Celsius temperature. Simultaneously, 
a Maxwell field TLM algorithm run in the same mesh provided the skin effect heat sources. 
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Fig. 2. Jet over backward facing step [ highly turbulent regime ]. 


Fig. 3. Karman vortex street behind cylinder[ snapshot of oscillating flow J]. 


T 

*mesh’ 
‘air_flow1? ——— 
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0.05 
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-0.1 
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——> 1.00000e-01 m/s 

-0.05 0 0.05 0.1 0.15 
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Fig. 4. Air convection in horizontal coaxial line RL100-230 ; 
transverse section [ reference arrow: 0.1 ms~'; transmitted power 160 kW CW, frequency 100 MHz; air dielec- 
tric; inner conductor copper, outer conductor cooled at 40 degrees Celsius ]. 
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Fig. 5. Heating process in coaxial line RL100-230 ; 
temperature maximum vs. time from power-on, 
convection compared to - fictitious - pure conduction 
[ transmitted power 160 kW CW, frequency 100 MHz, 
outer conductor cooled at 40 degrees Celsius J. 


Temperature profile in RL100-230 
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Fig. 6. Cross sectional temperature profile in coaxial line RL100-230 computed at steady state [ transmitted 
power 160 kW CW, frequency 100 MHz; air dielectric; inner conductor copper, outer conductor cooled at 40 
degrees Celsius ]. 
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